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We present an ab-imtio study of the phase transition cd— »/3-tin in Si and Ge under hydrostatic 
and non-hydrostatic pressure. For this purpose we have developed a new method to calculate the 
influence of non-hydrostatic pressure components not only on the transition pressure but also on 
the enthalpy barriers between the phases. We find good agreement with available experimental and 
other theoretical data. The calculations have been performed using the plane-wave pseudopotential 
approach to the density-functional theory within the local-density and the generalized-gradient 
approximation implemented in VASP. 

PACS numbers: 64.70.Kb 71.15.Nc 81.40.Vw 



I. INTRODUCTION 

The phase transitions in silicon (Si) and germanium 
(Ge) from the cd phase (cubic-diamond structure) 
to the /3-tin phase (body-centered tetragonal struc- 
ture, bet) are two of the most studied solid-solid 
phase transitions in condensed matter physics, both, 

experimentallyl234567891011121314151617181920212223 and 

thcorcticall; }'^^ ^^^^^^^"^' ^^^'^"'^^^^ ^^ '^'^'^'^"^'^^'^^'^^'^^^^"^""^^'^^ "^^ 
In the experiment, the phase transition in Si occurs at 
around 110 kbar and in Ge at around 105 kbar, where 
also lower values of the transition pressure are obtained. 
These lower values are often considered as caused by 
non-hydrostatic conditions, which are able to reduce the 
transition pressure, 

In fact, the pressure in the anvil cell is not ex- 
actly hydrostatic. Usually at pressures up to 100 kbar 
the pressure-transmitting medium yields nearly hydro- 
static conditions.'^- Above 150 kbar a non-hydrostatic 
pressure profile is visible, and at very high pressures 
the pressure-transmitting medium becomes solid which 
causes a strong non-hydrostatic effect. Even in the 
hydrostatic pressure regime there is a small pressure 
gradient 1^ Non- hydrostatic pressure profiles can also 
be an effect of the geometry of the celli^ Because of 
relaxation phenomena which happen in the pressure- 
transmitting medium the time for compressing and de- 
compressing has an influence on the measurement. 

In theoretical investigations generally hydrostatic con- 
ditions are assumed. Within calculations using the local- 
density approximation (LDA) the calculated transition 
pressures vary between 70 and 99 kbar for Si and be- 
tween 73 and 98 kbar for Ge. Usually the transition 
pressure is strongly underestimated by LDA calculations 
whereas calculations using generalized-gradient approx- 
imation (GGA) match the experimental value better 
(102-164 kbar for Si and 96-118 kbar for Ge). In any 
case, the discrepancy between experimental and theoret- 
ical results can also be due to non-hydrostatic pressure 
conditions in the experiment. Ab-initio calculations con- 
sidering non-hydrostatic pressure are rareSi54 and deal 
just with the transition pressure. The influence of non- 



hydrostatic conditions on the enthalpy barrier between 
the two phases is not studied within an ab-initio calcu- 
lation until now. Therefore, we developed a new method 
to calculate the transition pressure as well as the en- 
thalpy barriers also for non-hydrostatic conditions. In 
a first step, we obtain the transition pressure and the 
enthalpy barrier between both phases as a function of 
pressure starting from a complete numerical equation of 
:stff))i7e for hydrostatic conditions. Here a complete equa- 
tion of state means a continuous, multivalued function 
V{p) where V is the volume and p the pressure, simi- 
lar to the one of the common textbook example of the 
van der Waals gas. In a second step, this procedure is 
extended to non- hydrostatic conditions. 

This contribution is organized as follows: Firstly, a 
short introduction into the theoretical background of the 
calculations and the properties of the unit cell used for 
our calculations is given (Section ^ . Then we explain 
the procedure of calculating a complete equation of state 
from a given energy surface (Section llll|l . Following this, 
we present results for the influence of non-hydrostatic 
pressure components on the transition pressure and on 
the height of the enthalpy barrier (Section IIV|I . Finally, 
after a discussion of the results and comparison with 
available theoretical data (Section 0) we describe pos- 
sible extensions of our procedure and summarize (Sec- 
tion |Vlll- 

II. METHOD 

We have carried out calculations with the Vienna ab- 
initio simulation package (VASP)i^^*^SiSiS It is based 
one a plane-wave pseudopotential approach to the 
density- functional theory (DFT)i^24SS We have used ul- 
trasoft Vanderbilt-type pseudopotentials^ as supplied by 
Kresse and Hafner^^. The exchange-correlation poten- 
tial has been calculated within the GGA due to Perdew 
and WangSg for Si and Ge and the LDASlg^ for Si 
only. The forces on the atoms are derived from a gen- 
eralized formSiS& of the Hellmann-Feynman theoremSi 
including Pulay forcesiS^ For the ultrasoft pseudopo- 
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tentials a kinetic-energy cutoff of 270 eV (410 eV) for 
Si (Ge) has been sufficient for convergence of the to- 
tal energy and provides an error smaller than 0.5 kbar 
(0.2 kbar) for Si (Ge) to the pressure according to 
the Pulay stress:^ The special-point summation re- 
quired a 18x18x18 (24x24x24) mesh of Monkhorst- 
Pack points^^ which amounts to 864 (1962) k-points 
in the irreducible wedge of the Brillouin zone for Si 
(Ge). Since the /3-tin phase is metallic we have used a 
Methfessel-Paxton smearing2£ with a width of 0.2 eV, in- 
cluding the cd phase, since it is not a priori clear whether 
a given set of volume V and ratio c/a of lattice constants 
yields a metallic or a semiconducting phase. 

In order to minimize an energy offset between the 
structures it is important to describe the structures of 
both phases within the same bet cell (lattice constants 
a = b c) with two atoms in the basis at (0, 0, 0) and 
(0, 0.5a, 0.25c). The symmetry of the cd phase requires 
c/a = \/2, whereas c/a can vary for the /3-tin phase. Us- 
ing the bet cell the structure of the cd phase with respect 
to the conventional face-centered cubic cell is rotated by 
45° around the c-axis. In the following, CD and BCT 
denote the structure of the cd- and the /3-tin phase using 
the bet cell. 



III. COMPLETE EQUATION OF STATE 
UNDER HYDROSTATIC CONDITIONS 

Neglecting temperature and zero-point motion effects, 
the phase with the lowest enthalpy H = E + pV is sta- 
ble. Therefore, the transition pressure for a first-order 
pressure-induced phase transition from the cd phase to 
the /3-tin phase can be determined from the crossing of 
the corresponding enthalpy curves 

with H^'^{p*^) = ff^"*™(p*). First-order phase transitions 
are accompanied by a discontinuity in the volume AV 
and an overcoming of an enthalpy barrier which is lo- 
cated between the phases and which has a height of AH. 

Under hydrostatic conditions the pressure is defined as 
p — —dE/dV. It can also be determined from the stress 
tensor (Ti2i*2^ Since the structures here are orthogonal, 
the off-diagonal components of cr vanish, and cr has the 
form 




We are using a tetragonal cell, and thus Px = Py Un- 
der hydrostatic conditions all three components are equal 
and correspond to the external pressure p, 



Px^Py^Pz^P ■ (2) 

Under non-hydrostatic conditions the average pressure is 
defined as 

PO^\{Px+Py+Pz)^-\trCT , (3) 



which is again equal to the external pressure in the case 
of hydrostatic conditions. 




30 35 40 

Volume (A') 



FIG. 1: Contour plot of the total energy E{V,c/a) (solid 
lines) and of the average pressure po{V,c/a) (dashed lines, 
see Eq. for Si with GGA. The bold-dashed line corre- 
sponds to the value po — 0. The interval of the contour lines 
is 50 meV for the energy and 20 kbar for the pressure sur- 
faces. The black dots mark the equilibrium positions of the 
cd (c/a — \/2) and the /3-tin phase (c/a — 0.55). The dotted 
line marks the parameters under hydrostatic condition. 

We have calculated the total energy as a function of 
V and c/a. The corresponding energy surface E{V,c/a) 
is shown in Fig. ^ for Si using the GGA (similar re- 
sults are obtained for Ge within GGA and for Si within 
LDA). The two local energy minima, according to the two 
phases, with a saddle in between are visible. The pres- 
sure po{V, c/a) is obtained from Eq. and Eq. © and is 
included in the figure. Except along the (dotted) hydro- 
static line with px = Pz the pressure is non-hydrostatic. 
The local minima are at the crossing of the po = E^nd 
the hydrostatic p^ = Pz lines which defines the equilib- 
rium position. A similar graph can be drawn also for the 
enthalpy.^^ 

Along the hydrostatic line Px — Pz the structures are 
in a local equilibrium, meaning the total sum of internal 
and external forces caused by the pressure is zero. Hence, 
this condition can be used to extract the total energy E^ 
the external pressure po — Pi ^^nd the volume Y from 
Fig. ^in order to derive a full equation of state Vij)). 
Also iJ(p) and, furthermore, the values along the line 
across the saddle are accessible. These curves are shown 
in Fig. 13 Here, the ideal cd structure (c/a = \/2) has 
been reached only within an error of 1% for the lattice 
parameters. In order to discriminate the enthalpy curves 
against each other we have reduced them by a linear back- 
ground. The local stability is in accordance with the fact 
that the V(j)) curves are monotonically decreasing and 
the ff(p)-curves are convex. This is in contrast to the 
the textbook example, the van der Waals gas, where the 
line corresponding to the dotted line of i?(p) is concave 
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and signals local instability. 
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FIG. 2: Volume V, reduced enthalpy (see text), and enthalpy 
barrier AH as a function of the hydrostatic pressure for Si 
within GGA. The crossing of the solid and the dashed line 
determine the transition pressure. 

The transition pressure p* obtained from the cross- 
ing of the enthalpy curves are listed in Table U The 
corresponding change AV^ in the volume at the phase 
transition can be read from the upper panel of Fig. |21 as 
the difference between V'^'^{p^) (solid Hue) and T/'^~*'"(p*) 
(dashed line). Analogously the enthalpy barrier AH can 
be determined from the figures. In order to check the 
reliability of this method we compare the results with 
our previous ones^ based on the same total-energy cal- 
culations but obtained with a different method to evalu- 
ate the transitions pressures and enthalpy barriers. The 
agreement is very good, and the small differences are 
due to numerical errors. Thus, we can trust in the new 
method developed here. 



TABLE I: Transition pressures p*, volume changes AV, and 
enthalpy barriers AH derived from the complete equation of 
state in comparison with our previous results obtained with 
an alternative method (in parenthesis).'*'^ 





Ge-GGA 


Si-GGA 


Si-LDA 


p* (kbar) 


96 (96) 


122 (121) 


80 (79) 


AV (A^) 


7.5 (7.5) 


8.3 (8.3) 


8.5 (8.5) 


AH (meV) 


421 (423) 


510 (515) 


502 (508) 



we can calculate also the enthalpy barrier as a function 
of pressure. We have to distinguish between the bar- 
rier for the cd— >/9-tin transition, approaching from the 
cd phase, and the one for the /3-tin— >cd transition, ap- 
proaching from the /3-tin phase. In general, the enthalpy 
barrier AH has its origin in the energy saddle between 
the two phases and can be calculated as the difference 
of the reduced enthalpy of the phases (solid and dashed 
lines, respectively) and the one from the saddle (dotted 
line), see Fig.|21 In particular, the enthalpy barrier from 
the cd phase is the difference between the solid and the 
dotted line, and the one from the f3-tm phase is the dif- 
ference between the dashed and the dotted line. At the 
transition pressure p* the enthalpy barriers from both 
phases have the same height. The determination of the 
enthalpy barrier as a function of pressure is important to 
estimate the barrier in the case of over- and underpres- 
surizing the medium. Hence, the phase transitions will 
happen at a pressure different from p* which results in a 
different height of the barrier. As expected, the enthalpy 
barrier from the cd phase is decreasing with increasing 
pressure whereas the one from the /3-tin phase decreases 
with decreasing pressure. At zero pressure there is still 
an enthalpy barrier left for the /3-tin— >cd transition. This 
points at the fact that there is no spontaneous transition 
/3-tin— >cd. In the experiment the phase transition cd— >/?- 
tin is irreversible. 



IV. PHASE TRANSITION UNDER 
NON-HYDROSTATIC CONDITIONS 



The procedure for determining transition pressures and 
enthalpy barriers described in the previous section can be 
extended to the case of non- hydrostatic pressure. Besides 
the hydrostatic condition Pz — Px = the values for non- 
hydrostatic pressure components Pz — Px = d (with a 
fixed value of d) can be extracted from the energy surface 
E{V, c/a) along the corresponding lines of Fig. El A first 
naive trial considering just the total energy E^^^ under 
non-hydrostatic conditions and the corresponding values 
Pq^ for the average pressure and V^^ for the volume gives 
wrong results, e.g., an increase of the transition pressure 
for pz > Px and Pz < Px- This is in contrast to the ex- 
perimental observations. Thus is is necessary to include 
strain effects. 

Similar to the stress tensor of Eq. ^ the strain ten- 
sor e can be reduced to a diagonal form for orthogonal 
systems 



..71.72.73 



£22 



/ 



£33 / V 



(4) 



where Cx, £y, and are along the cartesian crystal axes. 
Since we have determined a complete equation of state, For small stress and homogeneous strain the components 
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FIG. 3: Contour plot of the total energy E{V, c/a) (solid 
lines) for Ge (GGA) with an interval of the contour lines 
of 0.2 eV. Besides the hydrostatic condition (bold solid line) 
non-hydrostatic conditions [pz — Px = — 15, —10, ... , 20) are 
shown. 



of e can be derived as^LS 



'^33 



(5) 



where Xj is the lattice parameter in the j-direction. Here 
Xj corresponds to the unstrained and x'- to the strained 
crystal. Including stress and strain the enthalpy can be 
written asS 



H = E 



(6) 



where H is the enthalpy and E the total energy per vol- 
ume. The calculation of the enthalpy at non- hydrostatic 
stress is based on Eq. The numerical realization is 
described in Appendix 1X1 

It turns out that the strain-only contribution to the 
enthalpy vanishes for the cd phase, that it is linear with 
the pressure for the /3-tin phase, and that it is non-linear 
for the contribution along the line across the saddle. This 
effect is apparent in Fig. 0] where the enthalpy including 
strain is presented. Since there is no strain effect on the 
cd phase, the change of the transition pressure is due to 
the strained (3-t'm phase. From Fig. 0] we can find the 
transition pressures for fixed non-hydrostatic conditions 
in the same manner as mentioned in the previous section. 

In addition to the average transition pressure pg also 
their components p* and in the x and z directions are 
shown as well as the enthalpy barrier at the phase transi- 
tion as a function of Pz — Px in Fig. [S] For the transition 
pressure one has the relation Pq — (2p^ +pl)/3. The 
boundary for the lowest pressure is fixed by the condi- 
tion that the components of po are not negative, since 
we exclude a stretching of the crystal. We find a strong 
lowering of the transition pressure if the pressure in the 
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FIG. 4: Equation of state -ff(po) for non-hydrostatic condi- 
tions as a function of the average pressure po. The difference 
Pz ~ Px of two neighboring lines is 5 kbar. The black dots 
mark the transition pressures pg. 
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FIG. 5: Enthalpy barriers at the average transition pressures 
Po and transition pressures (average transition pressure pg 
and the corresponding components p% and p' ) as a function 
of Pz — Px- The dashed line marks the boundary Pa; = and 
Pz = 0. 



z direction is larger than in the x and y directions. Thus, 
the crystal is more stable against compression along the 
X- and y-axes in the case of a strong non-hydrostatic 
component in these directions which cause an increase 
of the transition pressure. The corresponding enthalpy 
barriers are lowering in any case, but their value remains 
still larger than the thermal energy at room temperature 
(RT). 

Besides the non-hydrostatic effects we can consider fi- 
nally also the case of over- and under-pressurization of 
the crystal. To this end calculations of the enthalpy 
barriers as a function of the average pressure and non- 
hydrostatic conditions have been carried out. At very 
high pressures and very large non-hydrostatic compo- 
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FIG. 6: Enthalpy barriers like in Fig. |21 for fixed non- 
hydrostatic conditions Pz — Px = —20, —15, ... 35 kbar. The 
solid horizontal line at 25 meV marks the thermal excitation 
energy at RT. 



nents the enthalpy barrier for the cd^/3-tin transition 
is smaller than the thermal energy at RT, but these con- 
ditions do not appear by chance in the experiment; in- 
stead, they have to be applied by intention. In contrast, 
the enthalpy barrier is never smaller than 25 meV for the 
/3-tin— >cd transition even at the largest non-hydrostatic 
components with px, Py, and pz not negative (no stretch- 
ing). 



the linear functions which are = 0.651 -I- 35.1 kbar 
for Ge(GGA), = 0.606 p^-h47.3 kbar for Si(GGA), and 
pI = 0.619 -1-29.8 kbar for Si(LDA), respectively, agree 
very well with the ones of Refs. ^ and '54'. The differ- 
ence of the additive constants rely on the different values 
of the transition pressures in the hydrostatic case. Since 
Cheng et alm^^ restricted themselves to the enthalpy 
difference between the phases using path integrals, the 
enthalpy barrier was not accessible to them. 

The experimental values for the transition pressures 
vary between 103 and 133 kbar for SiL^i^il^ and be- 
tween 103 and 110 kbar for Go^i^i'^i^i^i^°i" where the 
firmed values are at around 110 kbar and 105 kbar, re- 
spectively. In both cases, our results obtained with GGA 
agree perfectly, whereas the LDA result underestimates 
the experimental value, which is a well known problem. 

The good agreement of our results with the ones of 
Cheng et alm^^ confirm the reliability of our method, 
which provides a larger field of applications. In ad- 
dition, our method can be extended to, e.g., the 
tin— >/mma— i'Sh transitions in Si and Ge. After the ex- 
traction of a two-dimensional energy surface from a three- 
dimensional one using the values along the lines where 
two components of the stress tensor are equal (like in our 
previous work^Si2&), the method mentioned here can be 
applied to this extracted surface. By the choice of two 
equal components the main pressure direction is chosen. 
Further extensions even to non-orthorhombic structures 
are possible, too. 



DISCUSSION OF THE RESULTS 



VI. SUMMARY 



In the past, non- hydrostatic conditions have been in- 
vestigated for different reasons i^2iSi2ii2Si2Si2I Directly 
comparable with our results are just the ones from Lee 
et alX'^ and Cheng et al^^^ Besides the transition pres- 
sures also the function p\{pl.) is given in these contribu- 
tions. This function can be obtained from our results 
by a linear fit of the values for the transition pressure 
in Fig. |S1 Already Lee et al7^ found a linear relation 
of p\ and p\ . In their molecular-dynamics investigation 
they found — 127 kbar for hydrostatic conditions and 
p\ — P\: ~^ 90 kbar. The additive constant corresponds 
to the lowest possible transition pressure. On the con- 
trary, Cheng et alM^M- obtained — 95 kbar and pi — 
0.737 pI + 25 kbar for Ge, and p* = 114 kbar and = 
0.658 pI. + 39 kbar for Si. Although the last values have 
been obtained also with VASP using GGA, the results 
for the transition pressure are slightly different from our 
results (95.9 kbar for Ge(GGA), 122.1 kbar for Si(GGA), 
and 79.6 kbar for Si(LDA)). This can be due to the fact 
that they have used different cells for the phases and also 
different pseudopotentials and convergence parameters. 
The choice of different unit cells can lead to an energy 
offset between the energy curves to which the transition 
pressure is very sensitive. Nevertheless, our results for 



We have developed a new method for investigating 
first-order high-pressure phase transitions which is based 
on the calculation of a complete equation of state. Be- 
sides the transition pressure and the volume change, 
which are also available with the common-tangent con- 
struction, the enthalpy barrier between the phases can 
be obtained with our method. A comparison with results 
for Si and Ge from common methods shows the reliabil- 
ity of the new method. Further on, the enthalpy barrier 
can be determined as a function of the external pres- 
sure which makes effects from over- and underpressuriz- 
ing accessible. An extension of this method allows us also 
to investigate high-pressure phase transitions under non- 
hydrostatic conditions, in particular the transition pres- 
sure and the enthalpy barrier, which are both decreasing 
if the pressure component along the c-axis is larger than 
the other ones. Our results show an excellent agreement 
with available experimental and theoretical data. This 
new method can be extended also to other phase transi- 
tions and also to ones including orthorhombic structures, 
for example, the transitions P-tin—^ Imma—f sh. Thus, we 
have developed a powerful tool for investigating phase 
transitions under hydrostatic and non-hydrostatic condi- 
tions. 
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APPENDIX A: CALCULATION OF THE 
ENTHALPY INCLUDING STRESS 

Here we give a short description of the formu- 
lae used for the calculation of the enthalpy includ- 
ing stress and strain effects starting from Eq. 

dH = dE 

( da dh dc\ (Al) 

+ \Px— +Py-r- +Pz— \Vo , 

V ao oo Co / 
which is equivalent to dH = dE + pdV in the case 



of hydrostatic conditions. Since Eq. lO holds just for 
small stress, we need to integrate this equation and can- 
not go directly to the absolute values. The integra- 
tion is performed using the recursively defined equation 

+ H'-^ 

where Xj are the lattice constants along the three carte- 
sian directions x^y and z, and the difference from the 
previous step (* — 1) is calculated along a line Px^Pz ~ d 
for fixed non-hydrostatic conditions (see Fig. O starting 
from the equilibrium structure of the cd phase. E^h is 
here the total energy along a line Px — Pz — d- The 
enthalpy H(p) under non- hydrostatic conditions corre- 
sponds to the calculated points H^(p^). By symmetrizing 
this equation numerical errors have been reduced. 
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